# Plot to illustrate regression line with confidence, prediction intervals
# 23 February 2013
# njg

default <- par(no.readonly=TRUE)
Ant.Data <- read.csv("antcountydata.csv",header=TRUE)
attach(Ant.Data)

##########################################
# fit the model

reg.model <- lm(mean.ann.temp~lat.centroid)

#################################################
# confidence and prediction intervals for original data

con.1 <- predict(reg.model,interval="confidence")
pred.1 <- predict(reg.model,interval="prediction")

####################################################
# confidence and prediction intervals for new x data
# (extrapolation and more finely spaced x)


nd <- seq(35,50,length=100)

con.2 <- predict(reg.model,interval="confidence",
                 newdata=data.frame(lat.centroid=nd))

pred.2 <- predict(reg.model,interval="prediction",
                  newdata=data.frame(lat.centroid=nd))

#####################################################
# plot results just for data interval

plot(lat.centroid,mean.ann.temp)
abline(reg.model)

matlines(Ant.Data$lat.centroid,con.1[,c("lwr","upr")],
         col="red",lty=1,type="l")

matlines(Ant.Data$lat.centroid,pred.1[,c("lwr","upr")],
         col="blue",lty=1,type="l")

##################################################
# plot results for extrapolated interval and smooth data
plot(lat.centroid,mean.ann.temp,xlim=c(30,50),ylim=c(-5,15))

abline(reg.model)

matlines(nd,con.2[,c("lwr","upr")],
         col="red",lty=1,type="l")

matlines(nd,pred.2[,c("lwr","upr")],
         col="blue",lty=1,type="l")

###########################################################
# plot the results as a confidence band using a polygon

plot(lat.centroid,mean.ann.temp,
     xlim=c(34,50),ylim=c(-5,25),type="n")

polygon(c(nd,rev(nd)),c(con.2[,2],rev(con.2[,3])),
        col="thistle",border=NA)

points(lat.centroid,mean.ann.temp)

abline(reg.model)

##########################################################
# plot the results as a prediction band using a polygon

plot(lat.centroid,mean.ann.temp,
     xlim=c(34,50),ylim=c(-5,25),type="n")

polygon(c(nd,rev(nd)),c(pred.2[,2],rev(pred.2[,3])),
        col="wheat",border=NA)

points(lat.centroid,mean.ann.temp)
abline(reg.model)

#########################################################
detach(Ant.Data)
par(default)
rm(list=ls())